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Abstract 

An analysis of adhesively bonded joints using conventional finite elements does not 
capture the singular behavior of the stress field in regions where two or three dissimilar materials 
form a junction with or without free edges. However, these regions are characteristic of the 
bonded joints and are prone to failure initiation. This study presents a method to capture the 
singular stress field arising from the geometric and material discontinuities in bonded 
composites. It is achieved by coupling the local (conventional) elements with global (special) 
elements whose interpolation functions are constructed from the asymptotic solution. 

Introduction 

Although bonded joints are a prime means for transferring load in the construction of 
composite structures, they are potential failure sites due to the presence of geometric and 
material discontinuities causing high stress concentrations. Reliable predictions of the gross 
response of the structure cannot be made accurately unless a precise description can be made of 
the interface through which the transfer of load is achieved. Thus, understanding the nature of 
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interfacial stresses is critical in designing reliable bonded joints, and efforts to understand the 
mechanisms needed to improve the strength of the bonded isotropic and composite materials are 
still continuing. Previous analyses of bonded joints can be categorized as the “shear-lag” and 
"finite-element” models. An extensive review and in-depth discussion of the previous 
investigations can be found in articles by Tsai and Morton (1994) and Ding and Kumosa (1994). 

Both the shear-lag models and finite-element models with conventional elements fail to 
capture the singular stress field at the junction of dissimilar materials. Blanchard and Watson 
(1986) concluded that a finite element analysis of such regions would not guarantee a convergent 
peak stress even with continued mesh refinement. In order to capture the exact nature of the 
stress field and to minimize the intensive computations arising from the refinement of the mesh, 
Barsoum (1988a, 1988b, 1990) introduced an iterative scheme in conjunction with the finite- 
element analysis without the use of a special (enriched) element. This approach is effective for a 
bimaterial interface with or without cracks. However, it suffers from the number of iterations 
required for convergence and the inability to enforce the continuity of traction components 
across the interface. Also, the rate of convergence and the accuracy of the results are dependent 
on the material properties and the scaling of the displacements during the iterations. Ding and 
Kumosa (1994) and Ding et al. (1994) applied this method to determine the singular stress field 
near the intersection of a bimaterial interface with free edges in adhesive joints. Although they 
captured the accurate description of the stress field near the junction, the strength of the 
singularity becomes inaccurate at distances very close to the free surface, where the failure 
usually initiates. This may be attributed to the limitation of the finite elements utilized in the 
analysis. 

To overcome this type of shortcoming in modeling a crack along a bimaterial interface, 
Chen (1985) developed an element with appropriate interpolation functions built in to account 
for the singularity at the crack tip. The unknown stress intensity factors are included explicitly in 
the expressions for the interpolation functions, and they are determined directly as part of the 
solution. Recently, Gadi et al. (1995) extended this method to determine the singular stress field 
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for the crack tip situated at the junction of three dissimilar sectors of material. Based on a 
similar concept, Kuo and Chen (1993) introduced a hybrid element with appropriate stress fields 
to investigate the transient thermal stresses in multi-layered regions with finite dimensions. 
However, both the hybrid and the enriched elements are limited to a specific geometry where the 
free edges are either perpendicular or parallel to the bimaterial interface, respectively. 

In a finite-element analysis of bonded joints, Destuynder et al. (1992) introduced a 
method to reduce the adhesive layer to a line through the use of asymptotic expansions of 
analytical solutions for the singular stress field near the geometric and material discontinuities. 
Also, Lin and Lin (1993) introduced a new element based on the Timoshenko beam theory for 
modeling the adhesive layer while accounting for the transverse shear and normal stresses in the 
adherents. However, this element does not account for the singular behavior of stress fields. 

The previous analytical and finite-element investigations were primarily concerned with 
isotropic adherent materials, because the presence of orthotropic materials is not suitable for 
directly constructing the analytical solution to the singular stress field. An analysis capability is 
lacking for determining the exact nature of the stress field in composite structures with bonded 
joints involving two or three dissimilar materials. Therefore, a global finite element with 
appropriate interpolation functions to capture the correct singular behavior arising from material 
and geometric discontinuities is developed here and is implemented into a finite element program 
with conventional elements. This program permits the analysis of various bonded joint 
configurations with accurate stress distributions in the critical regions and the extraction of the 
stress intensification parameter or the energy release rate in the presence of a crack. The results 
from this analysis provide the required parameters for the fracture criterion introduced by Gradin 
and Groth (1984) in order to estimate the strength of the bonded joint. Also, this program is 
integrated into the commercially available finite element program ANSYS so that the designer 
can use the ANSYS pre- and post-processing capabilities and execute the program within the 


ANSYS environment. 
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Solution Method 

The global-local finite element concept introduced by Mote (1971) is utilized in 
determining the stress field in regions consisting of a junction of two or three wedge-shaped 
sectors of orthotropic material with or without free edges (Figure 1). The extension and 
application of this method were demonstrated by Bradford et al. (1976, 1979, 1984), Dong 
(1983), and Her (1990). The development of the global element stiffness matrix is similar to that 
of a conventional (local) element, except for the interpolation functions. In this study, these 
functions are established by solving for the stress and displacement Fields in the regions 
illustrated in Figure 1. In these regions, each material is assumed to be elastic, homogeneous, 
and specially orthotropic, with elastic coefficient cjy . In reference to the Cartesian coordinates 
(x y y), the stress-strain relations under plane-strain assumptions are 



where <j ^ and are the stress and strain components, respectively. A method is presented 
in the Appendix that provides an average stiffness matrix for balanced laminates for which the 
material and reference coordinate systems do not coincide. This average stiffness matrix contains 
the independent coefficients of a specially orthotropic material. Throughout this study, the sub- 
or superscript k denotes a specific sector of the region. The interfaces among the adjacent 
materials, specified by angles 0, , are assumed to be perfectly bonded, thus requiring the 
continuity of traction and displacement components along the interfaces. 

The explicit forms of the displacement and stress components in the vicinity of the 
junction are constructed by solving for the displacement equilibrium equations expressed as 

C \\ u x.xx + C 66 Ll x,yy + ( C 12 + C 66 ) u y\xy = 0 
C 66 u y\xx + C [\ Ll y\yy +( C \2 + C 66 ) Ll x,xy =0 


( 2 ) 
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As suggested by Williams (1952), representing the displacement components as 

iia(r,8) = r*Ga(Q) with a = x,y (3) 

permits the reduction of Equation (2) to a system of ordinary differential equations in terms of 
the unknown functions Ga(8 ) in matrix form" 


M* ( 6 ■ C, )G l (0) + ( 1 - A)M* (0; C, )G' k (0) 


(4) 


A 2 M k (0; Cy ) + y (0; C tj )-X 1 G k {C ij ) G k (0) = 0 


where the known matrices and M k are defined by 


MCy) = 


(Cfi+c| 6 ) 


(c| 2 +c| 6 ) 


Mi = (0; Cij ) = 


m(0; , C*j ) 

—m'(9;Ci2,cj- 6 ) 


lm'(0;Cf 2 ,c| 6 ) 

m(0;C2 2 ,C6 6 ) 


with m{Q\ci,b) = a cos 2 0 + bsin 2 0. The unknown functions g£( 0) are contained in the vector 
G* ( 6 ) as G l ( 6 ) = [G* ( 6 ), G * (0)]. In Equation (2), the displacement components are defined 
in reference to a polar coordinate system, (r,0) , whose origin coincides with the junction of the 
vertices as shown in Figure 1. The unknown parameter A., dependent on the geometry and the 
material properties, indicates the strength of the singular behavior for the stress field. Utilizing 
the displacement representation given in Equation (3) and the strain displacement relations along 
with Equation 1, the displacement and stress components required for imposing the interface and 
boundary conditions can be expressed in polar coordinates as 


* Throughout this study, a prime denotes differentiation with respect to the variable 6 . 
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u k (r,d)=r X T(e)G k ( 6 ) 


<3 k(r, 6 ) = r X ~\XE k (0; Cy )G k (0) + F k (0; C i} )G' k (0)] 


in which the explicit forms of the matrices are given by 


T(0) = 


COS0 

sin0 


sin0 

COS0 


(5) 


m(d,C[ 2 ,C^ — 2Cgg)cos0 
-m(Q-,C \ i ~C \2 ~ C| 6 ,c| 6 )cos0 


—m(Q\ C22 + 2 Cgg , C\2 ) sin 0 

m (8’C66<C22 ~ C\2 -< ^66)sin 0 


and 


F*(0;Cy) = 


— m( 0 ; 2 + 2C56 , eft) sin 0 

m(0; Cge , C*] - Cy 2 ~ Cfojcosd 


m( 6 ,C 22 <Ci 2 + 2Cgg)cos 0 
-m(0; eft - C 22 - Cg 6 , Cg 6 , ) sin 0 


The vectors (r,0) and (r,0) contain the displacement and stress components, respectively, 

as u^(r,0) = [«* (r,0),«gj (r,0)] and c k =[cr QQ(r,6),G X g(r,6)] . For completeness, the normal 
stress component, a rr {r,6) is given by 

CT^(r,0) = r A_1 (Ae[ (0; C (j )G k (0) + f[ (0;Cy )G'* (0)] 

where 

e k (0; C,y) = [m(0; eft , eft +2C6g)cos0 m(0;Cft + 2 Cg 6 ,C 2 2 )sin 0 j 
f/T (0; Cy ) = [-m(0; C,*, - 2Cg 6 , eft ) sin 0 m(0; eft , C 2 * 2 - 2 C 6 * 6 ) cos 0 ] 

Determination of the unknown functions g£( 0) with a = x,y requires the solution of 
the system of ordinary differential equations with variable coefficients, Eq. (4), subject to the 
interface and boundary conditions. In the case of three sectors of dissimilar materials forming a 
junction, as shown in Figure la, the interface conditions are expressed as 
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°ap( r ’ 6 k) = V k ap( r ' d k) 

^a P (r,e k ) = <7 k ap 2 (r,e k+] ) 


k = 1,2 


> a,fi = r,9 with a = fi * r 


k = 3 


(6) 


“a( r - 0 iO = u a +1 (''’ 0 *) 


k = 1,2 


a,p = r,9 


ll a( r ’ e k)- u a l ( r ’9k+0 k =^ j 

For the intersection of a bimaterial interface with free edges, as shown in Figure lb, the interface 
and boundary conditions become 


<4/}O« 0 *) = 0 

o k a p{r,9 M ) = o k c ${r,9 M ) 
° k ap( r ’ e k + 2) = 0 


k = 1 


and a, [3 = r,9 with 


a = P* r 


(7) 


u a( r - 0 <:+l) = w a +1 (''' 0 A:+l) i k = ] and a,P = r,9 


The solution to the differential equations, Eq. (4), exists for values of \ that satisfy the 
characteristic equation of the homogeneous system of equations resulting from the imposition of 
the conditions (Eq. 6 or 7). Because of the complexity of the variable coefficients in Eq. (4), the 
solution to these equations is constructed numerically by recasting them as a set of first-order 
ordinary differential equations in terms of (0) and G'Jf ( 6 ) in the form 



G x (9) 

k 

0 

1 

0 

0 

k 

G x (9) 

d 

g; x (9) 


A| [ (0; Cjj , A) 

0 u (0;C y ,A) 

^12 fi< Gij , A) 

B 12 (0;C y ,A) 


G x (9) 

d9 

G y (9) 


0 

0 

0 

I 


Gy (9) 


G' y {9) 


A 2 \ (0; Cjj , A) 

0 21 (0;C y .A) 

<^22 Gjj , A) 

B 22 (9-C,j^)_ 


Gy{9) 


where A tJ and B tJ are the components of the matrices A and B defined by 


A^(0;C,,A) = M^(0;C,) 


A 2 M* (0; C, y ) + j M£ (0; C t] ) - A 2 L* (C y ) 
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and 

B* (0; C tJ , A) = ( A - 1)M* 1 1 (0; C i} )M' k (0; Q, ) 

The solution to this form of equations, Eq. (8), subject to the interface and boundary conditions 
is achieved by a Runge-Kutta forward integration scheme in conjunction with the shooting 
method. This procedure requires the initial estimates of the eigenvalues, X, and the corresponding 
eigenfunctions, G^(9 ) , as well as the target conditions. The integration process continues until 

the target conditions are satisfied, which requires the difference between the computed and 
prescribed conditions to vanish. In the case of three sectors of dissimilar materials, the target 
conditions are obtained from the interface conditions between the third and first regions as 

U 3 ( 03 ) = u,(0 4 ) and CT 3 (0 3 ) = o ] (0 4 ) (9) 

with 0 4 = 2k - 0 3 . The explicit forms of these conditions are expressed as 


Q = G 3 (0 3 )-G 1 (0 4 ) = O 


( 10 ) 

R = F 3 (0 3 ; Cy)G3(0 3 ) + A[E 3 (0 3 ; Cy) - E, (0 4 ; C,y )]G 3 (0 3 ) - F[ (0 4 )G{ (0 4 ) = 0 


in which Q T = [Q x ,Q y ] and R T =[R x ,R y ). In order to satisfy these conditions, Eq. (10), a 

positive definite objective function is defined in terms of the modulus of the complex functions 
Q x , Q y, R x , and R y as 

5=|e.r| 2 +|e v | 2 +|^| 2 +K| 2 (id 

With the well-established optimization techniques, the objective function is minimized by 
varying the real and imaginary parts of A, G^^^and G|( 04 ). 
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Determination of the eigenvalues and corresponding eigenfunctions satisfying the 
equilibrium equations and the interface conditions permits the expression of stress and 
displacement components as 


N 

G ap ~ 7 a p ( r > @ » A* i ) 

i= I 

N 

u a =Y, x ‘Ga( r ' e ^i) 

(=1 




( 12 ) 


The generalized coefficients, x f -, are determined by enforcing the continuity of the nodal 
displacements at the interface nodes between the global element and the surrounding local 
(conventional) elements. As illustrated in Figure 2, the global element with M interface nodes 
requires the imposition of the continuity conditions given by 


«x( r l.0|) 


g x {r\,Q\\X [ ) 

Gx( r \ ’0| ; Xi) 

1 

✓ — X 

qT 

C“ 

H 

' x l 

Uy(r { ,e [) 


Q y {r\,6\\X x ) 

Gy( r i ] . ^ i ; A 2 ) 

9y{r x ,6 x \k N ) 

x 2 

u y( r M '0 M ) 


Gx( r M ^ M'-X [) 

Gx( r M '9 M'Xi) ' 

' Gx( r M M 


u y (r M ,9 M ) 


G y( r \l 'Q M'<X\) 

Gx( r M M 'X 2 ) • 

" G y( r M M'X m\ 

X N 


or 

{«}=[$]{*} 

In general, the number of equations, 2Af, exceeds the number of unknown coefficients, A, 
resulting in an overdetermined system. Therefore, the unknown coefficients are expressed in 
terms of nodal displacements based on the least squares minimization procedure as follows: 

{.*}=[Z]{«} with [Z\ = [[Q] T [Q]]' [Q] T 


(14) 
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Determination of these coefficients permits the expression of the stress and displacement 
components in terms of the nodal displacements, and the strain energy in the global element 
becomes” 

U=±j[u} T [Z] T {7 al5 ) T {g a }[ZUu}ri j3 dS, a,p = x,y (15) 

5 


where Tj p are the components of the unit normal to the surface, S, of the global element. The 
vectors [J a p } and [§ a } are defined as 


{%}={?«/? ('-Ml), 7 a p(r,e-,X 2 ), .... J a/3 (r,0; A, v )} 


(16) 


Minimizing the strain energy with respect to the nodal displacements associated with the 
global element results in the global stiffness matrix [k] defined as 


[k] = ^j[Z}{T a p} T {g a } + {g a } T {T a p})[Z}TlpdS (17) 

5 

The global and local element stiffness matrices are assembled to establish the system equilibrium 
equations as 

{K]{5} = {F} (18) 

where [K] is the system stiffness matrix and the vectors {5} and {F} include the total nodal 
displacement and force components, respectively. This process led to the development of a finite 
element program incorporating both global and local elements. The local elements consist of 
quadrilaterals and triangular elements, whose interpolation functions can be found in any 
textbook on the elementary finite element method. The shape of the global element can be an n- 
sided polygon, depending on the details of the mesh surrounding the global element. The 
number of interface nodes and the size of the global element were established based on 
convergence requirements. 

” Repeated subscripts imply summation. 



This global-local finite element program is implemented in the ANSYS platform through 
the use of ANSYS Parametric Design Language (APDL) commands. It permits execution of the 
program without leaving the ANSYS environment. The global element is introduced into the 
ANSYS element library as “USER 104” through the User Programmer Features (UPFs) routines. 
For this purpose, ANSYS is customized and relinked by including two FORTRAN routines, 
“uecl04.f ” and “uell04.P\ The first routine describes the characteristics of the element, such as 
the maximum number of nodes associated with the global element. The second routine arranges 
the element matrices, load vectors, results, and element solutions data during normal ANSYS 
execution. However, it acts as a dummy routine because the global-local analysis program is 
used in constructing the solution rather than ANSYS. The real constants for each global element 
are defined by: 

1. The x and y coordinates of the origin of the local coordinate system associated with 
the global element. 

2. Material number associated with the region. 

3. Angle specifying each region. 

This capability permits the use of ANSYS pre- and post-processing for the global-local finite 
element analysis. 


Numerical Results 

Analysis of bonded dissimilar composite materials by the present approach is 
demonstrated through a single-lap adhesive joint with three typical configurations. The geometry 
and dimensions of each of these lap-joint configurations are described in Figure 3. The 
parameters C\ and c 2 and /zj and /z 2 denote the end distance (c) and thickness (h) of the top and 
bottom adherents, respectively, with numerical values of c\ = c 2 = 200 mm and 
h\ = h 2 = 5 mm . Also, a long and a short joint with overlap lengths of l = 320 mm and 40 mm 

are considered in order to capture the effect of joint length. For joint type II, bevel angles of 
0jand 0 2 are equal and are specified as 45°. The parameters f 2 , 0i» and 0 2 describing the 
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overflow of the adhesive in joint type m are specified as = i 2 = 2/i,with the adhesive 
thickness h = 0.4 mm and <p j = 02 = 45°. As shown in Figure 3, the upper adherent is subjected 
to a uniform stress, o*q, and the lower adherent is fixed at the other end. The adhesive is an 

isotropic material with Young’s modulus E = 3400 MPa and Poisson’s ratio v = 0.35. Top and 
bottom adherents are composed of [0°/90°/0°] plies with properties E L = 147 GPa, 
Ej =llGPa, G u =5.3GPa, and v u =0.3. The averaged orthotropic properties for the top 
and bottom adherents are computed to be 1 18.73 = GPa, C 2 2 =36.14 GPa, C [2 = 12.06 

GPa, and C 66 =6.21 GPa. In the case of isotropic adherents, the Young’s modulus and 
Poisson’s ratio are taken as E = 200 Gpa and v = 0.3, respectively. 

The finite element representation of each lap-joint configuration with global and 
conventional elements is illustrated in Figures 4-6. Their overall deformations under the 
specified load of 0*0 = 1 MPa are shown in Figure 7. The eigenvalues retained in the 

construction of the interpolation functions for each global element are tabulated in Table 1 for 
the isotropic adherents. For a type I joint, the behavior of the peel and the shear stresses in each 
global element along the bond line from the junction point is given in Figure 8. The peel and 
shear stresses along the bond line of the most stressed region represented by global element D for 
all joint types are shown in Figure 9. 

Based on the strain energy density criterion introduced by Sih and Macdonald (1974), an 
examination of the strain energy density around the junction point for a specified distance 
provides possible failure sites and the crack propagation path once the failure initiates. The 
variations of the tangential and shear stresses and the strain energy density for a specified core 
region, =0.05 mm, around the junction in each global element are shown in Figures 10-12. 
These figures reveal that the failure in joint type I is most likely to initiate in global element D at 
the junction of the top adherent and the adhesive, and it is predicted that it will propagate along 
the bond line. For type II lap-joints, the possible failure site is also in global element D along the 
bond line. For type III lap joints, the failure may initiate at the junction points E or F and the 
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crack is likely to grow into the adhesive in the vertical direction. These results are all based on 
the eigenvalues presented in Table 1. 


Conclusions 

The global-local finite element analysis eliminates the use of a fine mesh and provides an 
accurate description of the stress field in the critical regions of the bonded joints. The order of 
the singularity along with the corresponding stress intensification parameter can be used for 
predicting failure in the adhesive layer of the joint. With this capability, the geometry and 
material properties can be optimized to minimize stress intensification. Also, this approach can 
be extended to the analysis of bonded joints with visco-elastic adhesive layers. As indicated by 
Ratwani et al. (1982), the effect of modulus relaxation becomes important because a large 
redistribution of stresses occurs while the joint is loaded. 
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Appendix 

Because only specially orthotropic materials are considered in this formulation, it is 
limited to laminates with a ply orientation of either 0° or 90°. However, a method exists to 
determine the average specially orthotropic stiffness matrix for a balanced laminate. A balanced 
laminate is a panel that has a negatively oriented ply for every positively oriented ply. Thus, this 
analysis can be expanded to model any balanced laminate if an average stiffness matrix is 


utilized. 
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To determine the average stiffness matrix, the stiffness matrix of each ply in the global 
coordinates is required. The stress-strain relationship for the k* layer of a laminate can be 
represented as 

o k = C k z k (19) 


or 

5 k =C k z k (20) 

where the tilde denotes the quantities in the local reference frame. The unit vectors of the local 
reference frame can be written in terms of the direction cosines and the global unit vectors as 



k 

~h 

m\ 

"l 


n x 

n n 

> = 

h 

m 2 

n 2 


n V 

A. 


h 

m 3 

n 3 . 


n. 


( 21 ) 


With these direction cosines, stresses and strains can be transformed to the global reference 
frame through the transformation 


_ k rp k k 

0=10 

t k = T k e k 


(22) 


where 
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(m 2 n 3 + m 3 n 2 ) 

( m l n 3 + m 3 n 0 

{myn 2 + ^ 2 ^ 1 } ) 

l\ n \ 

l 2 n 2 

h n 3 

(h n 3 + h n 2 ) 

(h n 3 +l 3 n l ) 

(/,n 2 +/ 2 «i) 

lym 

l 2 m 2 

l 3 m 3 

(l 2 m 3 + / 3 m 2 ) 

(/,m 3 +/ 3 m,) 

(/l m 2 +/ 2 " 1 !) 


Substituting Eqs. (21) and ( 22 ) into (19) while noting that 

'j'k h j 

yields 

a k =T k C k T kT e k 


(24) 


(25) 
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Thus, the transformation of the stiffness matrix from local coordinates to global coordinates is 

C k =T k C k T kT (26) 

The average stiffness matrix, C, of a balanced laminate is then represented by 

C=4iVc* (27) 

h 1 

which represents an average weighting of the stiffness matrix based on the thickness of the layer. 
Provided the laminate is balanced, C will represent a specially orthotropic material with 12 
nonzero coefficients, 9 of which are independent. This averaging process is expected to provide 
reasonably acceptable results, provided the paired balancing plies are located closely to one 
another and the plies are relatively thin. 
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Table 1. Eigenvalues associated with the global elements in the joints. 
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Figure Captions 

1. Junction of three and two dissimilar orthotropic materials. 

2. Interface nodes between the global element and the surrounding local elements with 
or without free edges. 

3. The geometry of a single lap joint: (a) without adhesive overflow; (b) with beveled 
adherents; (c) with adhesive overflow. 

4. Finite element discretization of the lap joint without adhesive overflow. 

5. Finite element discretization of the lap joint with beveled adherents. 

6. Finite element discretization of the lap joint with adhesive overflow. 

7. Overall deformation of the long and short joints. 

8. Variation of stresses along the interface near the junctions of a short and long joint 
without adhesive overflow in global elements A-D. 

9. Variation of stresses along the interface in global element D near the junction of 
short joints of type I-III. 

10. Variations of the tangential and shear stresses and the strain energy density around 
the junctions of short joints without adhesive overflow. 

12. Variations of the tangential and shear stresses and the strain energy density around 
the junctions of short joints with beveled adherents. 

13. Variations of the tangential and shear stresses and the strain energy density around 
the junctions of short joints with adhesive overflow. 
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